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Abstract 

pH . We investigate a high-order homogenization (HOH) algorithm for periodic mul- 

tilayered stacks. The mathematical tool of choice is a transfer matrix method. 
Expressions for effective permeability, permittivity and magnetoelectric coupling 
are explored by frequency power expansions. On the physical side, this high- 
order homogenization uncovers a magnetoelectric coupling effect (odd order ap- 
proximation) and artificial magnetism (even order approximation) in moderate 
contrast photonic crystals. Comparing the effective parameters' expressions of a 
stack with three layers against that of a stack with two layers, we note that the 
magnetoelectric coupling effect vanishes while the artificial magnetism can still 
be achieved in a center symmetric periodic structure. Furthermore, we numeri- 
cally check the effective parameters through the dispersion law and transmission 
curves of a stack with two dielectric layers against that of an effective bianiso- 
tropic medium : they present a good agreement in the low frequency (acoustic) 
band until the first stop band, where the analyticity of the logarithm function 
r — . of transfer matrix (log{T}) breaks down. 

oo 

^". ' 1. Introduction 

^^ ■ There is a vast amount of literature in the homogenization of periodic struc- 

tures with the classical effect of artificial anisotropy. A less well-known effect is 
artificial magnetism and low frequency stop bands through averaging processes 
in high-contrast periodic structures [f|, |2|, y, |4| ■ 

In layman terms, O'Brien and Pendry observed in 2002 that rather than 
using the LC-resonance of conducting split ring resonators to achieve a nega- 
tive permeability as an average quantity [5| , one can design periodic structures 
displaying some Mie resonance, e^g. by considering a cubic array of dielectric 
spheres of high-refractive index [lj. These resonances give rise to the heavy- 
photon bands in photonic crystals which are responsible for a low frequency 
stop band corresponding to a range of frequencies wherein the effective per- 
meability takes extreme values [3j. So-called high-contrast homogenization [4J 
predicts this effect which is associated with the lack of a lower bound for a 



frequency dependent effective parameter deduced from a spectral problem re- 
miniscent of Hclmholtz resonators in mechanics. 

In this paper, we would like to achieve such a magnetic activity without a 
high-contrast material. The route we propose is based upon a homogenization 
approach for high-frequencies i.e. when the period of a multilayered structure 
approaches the wavelength of optical wave. The extension of classical homogeni- 
zation theory [y, |7|, |8j to high frequencies is of pressing importance for physicists 
working in the emerging field of metamaterials, but applied mathematicians also 
show a keen interest in this topic @, 9, H|, where the periodic structures at 
sub-wavelength scales (A/10 to A/6) l|, [5| can clearly be regarded as almost 
homogeneous. 

The tool of choice for our one-dimensional model is the transfer matrix me- 
thod, which allows for analytical formulae as shown in Section [2] a high-order 
homogenization (HOH) method is proposed wherein Baker-Campbell-Hausdorff 
formula (BCH, an extension of Sophus Lie theorem) was implemented ; we stress 
that ideas contained therein can be extended to two and three-dimensional per- 
iodic structures, at least for simple geo metries, such as woodpile structures, of 
particular interest in photonics, see [11(. Importantly, we not only achieve ma- 
gnetic activity in moderate contrast dielectric structures, but also unveil some 
artificial bianisotropy in Section [3J where a multilayered stack with an alterna- 
tion of two dielectric layers is considered. As proposed by Pendry, the bianiso- 
tropy is yet another route towards negative refraction |12j, however, it is noted 
that the artificial bianisotropy vanishes in a periodic stack with center symme- 
try, where an extension of the HOH method applied to a stack with m(> 3) 
alternative layers is explored in Section [4] Furthermore, a correction factor is 
investigated in Section[5]to estimate the asymptotic error, which is proportional 
to l/n p with p the approximation order. We then numerically check the effec- 
tive parameters through the dispersion law and transmission curves between 
the multilayers and the effective medium in Section [6] : A good agreement in 
the low frequency band up to the first stop band verifies the equivalence of the 
two structures. We finally take a frequency power expansion of the transfer ma- 
trix, which is analytic in the whole complex plane in Section and draw some 
conclusions in Section [5] 

2. Mathematical setup of the problem 

Figure (Ha) shows a schematic diagram of periodic multilayered stack with 
a unit cell made of two homogeneous dielectric layers C\ and £2 of respective 
thicknesses h\ and hi (h = hi + hi = hi/n + h%/n), where n is the number of 
the unit cells in the stack. The permittivities of the two layers are assumed as 
Ei, £2 , and the permeabilities are /ii = fJ,2 = Mo- It should be noted that the 
whole thickness of the stack h(= nh) is a constant and comparable with the 
wavelength of light as h/X ss 1. Starting from the case when n = 1, the stack 
is a most simple structure consisting of only two dielectric layers, the effective 



medium theory [13j, l!4|, l!5( cannot be applied when h/X ~ 1. However, if we 



increase n to a large enough constant, then the system contains n times smaller 



unit cells, the thickness of which will be much smaller than the wavelength 
of light (e.g. h/X <C 1). Hence a homogeneous medium with permittivity e e g, 
permeability [i c g and bianisotropy K c g shown in figure [TJb) can be assumed 
to behave as an effective medium for such a multilayer. The approximation of 
the multilayer by the effective medium will be more accurate for large n, a fact 
which will be proved in the following sections. 
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Figure 1: (a) Schematic diagram of a multilayered stack consisting of an alternation of two 
homogeneous dielectric layers of permittivities ei, £2 and thicknesses h± = hi/n, hi = h\/n, 
with n the total number of the unit cells. The whole thickness of the stack is denoted by 
h = nh = n(hi + hi), which is comparable with wavelength by h/X ss 1, when n is large 
enough, then h/(n\) = h/X <C 1, in other words, uih/(2nco) <C 1 with co the velocity of light 
in vacuum, (b) An effective medium described by anisotropic tensors of permittivity, permea- 
bility and bianisotropy (i.e. a metamaterial with artificial bianisotropy and magnetism), the 
thickness is h. 



2.1. Time-harmonic Maxwell's equations 

At the oscillating frequency u, the electric and magnetic fields E and H 
are related to the electric and magnetic inductions D and B through the time- 
harmonic Maxwell's equations, 

V x H(x) = -iwD(x) , V x E(x) = iwB(x) , (1) 

and the constitutive relations for non-magnetic isotropic dielectric media, 

D(x) = e m E(x) , B(x) = H(x) , x 6 C m , (2) 

where e m the permittivity in the m homogeneous layer located in the domain 
C m of M 3 . 

Then a Fourier decomposition is introduced for both electric and magnetic 
fields as 



E(fei, k 2 , x 3 ) = — / ~E(xi,x 2 ,x 3 ) exp [ - i(kiXi + k 2 x 2 )] dx\dx 2 , 
H(fci, k 2 , xa) = — / H(x 1 ,x 2 ,x 3 ) exp [ - i(k 1 x 1 + k 2 x 2 )} dx 1 dx 2 . 

^ ./B2 



(3) 



with fci, &2 the projections of wave vector k on X\, x 2 axes, respectively, where 
an oblique polarizable plane wave with wave vector k = kixi + k 2 x 2 + ^3X3 is 
considered {xj, axis is perpendicular to the layers). 

Applying the decomposition © to equations ([T]) and ©, we derive an ordi- 
nary differential equation (involving 4 x 4-matrices and a 4-components column 
vector) [16J 



OF 

- — (w,fci,fc 2 ,z 3 ) = iM m (iu,k l7 k 2 )F(uj,k 1 ,k2,X3) 
0x3 



(1) 



where F is a column vector containing the tangential components of the Fourier- 
transformed electromagnetic field (E.H), i.e. the components along x\ and x 2 
axes. 

Here, in order to apply a much simpler notation for subsequent calculation, 
we would like to define a new set of coordinates denoted by (x\\ , x±, X3) : where 
the component x\\ is along the direction of wave vector k = (/ci, k 2 ), x± is along 
k = {~k 2 ,ki) which is perpendicular to Xn. In other words, the new set of 
coordinates is a rotation of the previous coordinates around the X3 axis. For 
every vector x, the change of coordinates from (x\, x 2 , x 3 ) to (xii, x±, x%) can 
be expressed by 
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Note that, thanks to the symmetry of the geometry, the parameters of the 
multilayer are invariant under this transformation 17 1. 

Hence, (j4]) can be recast in the new coordinate system as 



dF 
dx 3 

with the column vectors 

F11 



(u,k,xs) = iM m (u,k)F(u,k,x 3 ), 



F = 



F„ = 



E„ 



Hi, 



F, = 



Ei 
H 



(6) 



(7) 



where E||, H|| (respectively Ej_, Hj_) are the components of the electric and 
magnetic fields along the x\\ axis (respectively the x± axis). 

Correspondingly, the matrix M m is a 4 by 4 matrix, the components of which 
can be expressed as 



M m (u,k) =uj 







a m + {k 2 /u 2 )a m l 




/'n 




(8) 



where k = uk is the projection of wave vector k on atrii , and k — svnOi, with 
the incident angle. 



The matrix M m (u>, k) is independent of x 3 in each homogeneous layer, the 
solution of the equation ([6J in the layer C m is simply 



F(w, fc, £3 + h m ) = exp [iM rn (uj,k)h m ] -F(uj,k,x 3 ) 



(9) 



The exponential above is well-defined as a power series of matrix M m (u, k), and 
defines the transfer matrix in the m th layer of thickness h rn . Since this power 
series has infinite radius of convergence, the transfer matrix 



T m (uj, k) = exp[iM m /i TO ] 



(10) 



is analytic with respect to the three independent variables u), k\ and /c2- For an 
arbitrary permittivity profile (with the classical assumption of upper and lower 
bounds of permittivity greater than e uniformly in position x and number of 
layers n), analyticity is proved using a Dyson expansion [18j |. 

2.2. Main homogenization result 

From physical considerations, for the effective medium in figureQJb), we pos- 
tulate that the homogenized constitutive equations emerging from the asympto- 
tic limit n — > +00 in the sequence of equations ([2]) in the new coordinate system 
are 

D eff (fc,x 3 ) = e e ff(w, fe)E eff (fe, #3) +iK e s(u>,k) JH e e(k,x 3 ) 

B off (fc, x 3 ) = /U e ff(w, fc)H cff (fe, x 3 ) + iJK cS (uj, fc)E cff (fc, x 3 ) 

where s e s, /i g are tensors of rank two which represent respectively the (aniso- 
tropic) effective permittivity, permeability 
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matrix J corresponds to the 90 degrees rotation around the x 3 axis, and K c ff is 
the bianisotropic parameter measuring the magnetoelectric coupling effect 
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(13) 



Now, applying equation ([3]) to (Q]) and (|11|) . we obtain 



with matrix 



dF 

- — (u,k,X3) = iM cS (uj,k)F(uj,k,x 3 ) , 

dx 3 



M eS (u), k)=u 



-icr' K <t± + (k /uj 2 )^ 
-U|| ia K 



(14) 



(15) 



where the 2x2 blocs are defined by 
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(16) 



(17) 



These parameters are all unknowns at this stage which we would like to derive 
from a homogenization algorithm. The transfer matrix is correspondingly, 



T e ff(w, k) = exp[iM cB h\. 



(18) 



3. High-order homogenization (HOH) algorithm for multilayered stack 

Since we have derived the transfer matrix for the multilayered stack and 
postulated its structure for the effective medium in the previous section, it 
follows that the description of the homogenization procedure shown in figure [1] 
can be expressed as 



exp[iAf2^2] exp[iMi/ii] = ex-p[iM e gh] . 



(19) 



This means the two structures should present a same transmission property, 
where M e g is the unknown to be calculated. Note that the left side of equation 
(fT9|) is a product of two exponential functions, which can be approximated 
by introducing the Baker-Campbell-Hausdorff (BCH) formula (an extension of 
Sophus Lie theorem, see (l9|). In mathematics, the BCH formula is concerned 
with 

exp[Z] = exp [At] exp [A 2 ] (20) 

with A\ and A 2 square matrices. An analytical expression for Z is 



Z = log(exp[,4i]expL4 2 ]) 
1, 



A ± + A 2 



Mi,A a 



12 



{A^lAi.M 



12 



\A 2 ,\A U A 3 



(21) 



where \A\ , A 2 \ — A\A 2 — A 2 A\ is the commutator of A\ and A 2 , the product of 
which is noncommutative with [A2,^4i] = — [Ai, A2]. Here, we would like to de- 
note Ai+A 2 in (|21|) as the zeroth order approximation for Z, which corresponds 
to the classical homogenization ; |Ai, A2J/2 the first order, |^4i, [Ai, A2H/I2 — 
\A 2 , [Ai, A2U/I2 the second order approximation, and so on. 
From (|1T|) and (fl9l. we have 



Mi/ii 



1 



iM cS h = i(M 2 h 2 ,..,.,. 

1 
+ —fiM 2 h 2 , liM 2 h,2,iM 1 hi 



12 l 



[*M 2 /i2,iMifti] 

^{iMihuliMahaJMihi] 



(22) 



Furthermore, the expressions for the effective parameters in (fT2|) and (fill)) can 
be derived by comparing the two matrices in the left and right hand sides of 

(E2D- 

First, we consider the zeroth order approximation in (|22p . it yields M e g « 
M1/1 + M2/1 with the filling fractions fi(— hi/h) and f 2 {= h 2 /h), respectively, 
and the effective parameters are 



(23) 



£|| =e± = £1/1 + £2/2 , £3 1 = £1 Vi + e 2 V2 
A»|| = Ml = M3 = Mo , #|| = .Kj. = . 

They are identical to the effective permittivities presented in [1J, |l5|, |20|, |2l|, [22J 
by classical homogenization : the effective permittivity, permeability are equal 
to the average of two dielectric layers, while the bianisotropy is zero. 
If we go further by taking the first order approximation, we obtain 



M cS « A'hh + M 2 f 2 + y [M 2 / 2 , M1/1] 



(24) 



Since both Mi and M2 are off-diagonal matrices, then their commutator leads 
to a diagonal matrix, the components of which correspond to those of M c g in 
(USD, i.e. 



y[M 2 / 2 ,Afi/i| =w 



'K 



where 



, £1 - £2 , , 

CTK = W/l /1/2 
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k £i+£ 2 
W 2 £i£ 2 



(25) 



(26) 



provides the first order correction to the leading order approximation (classi- 
cal homogenization) in (|23p . This first order correction is encompassed in the 
following bianisotropic parameter 



K ± (w, k) = — Mo(ei - £2)/i/2 

.2 



#11 (<*>>&) = -r-Mo(ei-ea)/i/2 



1 



fe £1 + £2 



U^ /i £i£2 



(27) 



Notice that if is not only frequency dependent but also exhibits spatial disper- 
sion. It leads to if 1 1 7^ K± when k ^ 0. 

Furthermore, if we consider the second order correction, a term with "double 
commutator" will appear in the asymptotic expansion 



ih 



M cS » (M 2 / 2 + M1/1 

+ ^[M 1 / 1 ,[M 2 / 2 ,M 1 / 1 



[M a /a,Mi/i 



12 



IMa/a.IMa/^Mi/! 



(28) 



with commutator of Mi and M 2 given in (|25[) , and double commutator 



IMi/i,[M a / 2 ,Mi/i 



2w 2 



-/i 



fc 

(TiiTr- + CT^-CTi + — „- (Oi (7^ + v' K V\ ) 

o-i a' K + a K a x 



According to the definitions of a m in (|BJ| , or- and ct^- in (|1T[) , we have 
Thus, (I2U1) can be simplified to 



4w 2 
[Mi/i,[M 2 / 2 ,Mi/i]] = — /i 



fe -i 
uicrfiT H — t^ ctr: 



Or (Tl 







(29) 
(30) 

(31) 



Similar equalities hold for [M 2 / 2 , [Af 2 / 2 , A/1/1]]. Hence, the terms arising from 
"double commutator" lead to 



12 



(IMxA, [M 2 / 2 , Mi/!]] - [M 2 / 2 , [Ma/a.Mi/x]]) 



L0 2 h 





CK&lfl ~ CTk(T 2 / 2 



cifif/i - cr 2 <7x/ 2 H -(o^ (Ta'/i - ovT ctr-/ 2 ) 





(32) 

which is again an off-diagonal matrix. Substituting (|3"2"|) into (|28p and comparing 
with the form of matrix M e g in (|1 5|) , we find 



T \\ 



&ifi - 02/2 + -z-i&K&ifi — cr K a 2 f2) 



C-L = CT1/1 + (7 2 / 2 + —{cri(TKfl - <T2<TKh) 

Furthermore, the expressions of effective e e g, /x e ff are as follows 

, .2l2 / r>2 p 1 _ 

£||(w,fe) =£1/1 +£2/2 + — — IM»hfc( £ i -£2)(ei/i -£2/2) ( 1- 7-2— — 



(33) 



ex(w,fe) = £1/1 + £ 2 / 2 
ui 2 h 2 
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or ,uo£ie 2 
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M||(w,fe) = Mo - 
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(34) 



All the effective parameters are frequency dependent and with spatial dispersion. 
These expressions turn out to be equivalent to the ones reported in [23J, [2l|, |24[ , 
where the effective refractive index n e s is expanded using a power series of 
period-to- wavelength ratio A/ A. Taking equation (5)(s-polarized incidence is 
considered) in paper of Yeh [24] as an example, the dispersion relation for a 
two-component layered medium is approximated by taking the fourth order of 
0[(A/A) 2 ] as 

K* +0* * {^-f + {ab)2{nl ; 2 f HUj/c)4 05) 

where K and /3 are the z and x components of the Bloch wave vector, a and b 
are the thicknesses of the alternating layers, A = a + b is the period, n\ and 112 
are the indices of refraction of the corresponding layers, c the velocity of light 
in vacuum, and 



hi, b= h 2 , A = h, n 2 = — n\ 



(36) 



(37) 



2 2 2 

n o = J n l + J n 2 

Comparing with the notations in our formula, we have 

K = k 3 , (3 = k, n\ = e\, n\ = e 2 

o 

A 7 ^ 1 ' A' 
hence /i = a/ A, fa — b/A, then ([35]) is 

2 4 

kl + fc 2 « ^( ei A + e 2 / 2 ) + j^/i 2 /l(ei - £2) 2 (38) 

which contains the terms of w 2 and a; 4 . On the other hand, for the effective 
medium in our HOH process, the dispersion relation of k% versus uj is 

2 

k\ = — (ej_M|| - Kj_) -k 2 , s - polarization (39) 



Let us substitute the expressions of effective parameters (|34|) into (|39[) . and 
collect the terms up to w 4 in the calculation process, and finally we obtain the 
same formula as in (|38|) . Similar calculation applied to a p-polarized incident 
wave shows again our homogenization provides exactly the same effective index. 
Hence, it is stressed that the effective parameters e c ff> Meff a- n d K e g achieved 
from the HOH algorithm contain more information than the single effective 
index parameter in [2lL |24| . [23], e.g. the artificial magnetism and bianisotropy 
from periodic dielectrics which can not be seen in the derivation of refraction 
index. 

In the asymptotic process, we have noticed that the magnetoelectric coupling 
comes from the odd order approximation while the artificial magnetism and high 
order corrections to permittivity emerge from the even order approximation in 
([2i"j) . This can be explained in the following way : The matrix M m (m = 1,2) for 
the dielectric layer is off-diagonal, the terms of odd order approximation usually 
contain odd commutators, hence a diagonal matrix will result, the components of 



which correspond to K. However, the terms of even order approximation contain 
even commutators, the resulting matrix is always off-diagonal. This introduces 
the artificial magnetism and high order corrections to permittivity. 

Moreover, these results are fully consistent with descriptions in terms of 
spatial dispersion [25|, [26| where, expanding the permittivity in power series of 
the wave vector, first order yields optical activity and second order magnetic 
response. The equivalence of these two descriptions (frequency and wave vector 
power series) is confirmed by considering a unit cell with a center of symmetry, 
for example a stack of three homogeneous layers (permittivity e m and thickness 
h m , m = 1, 2, 3) with e 3 = e\ and h 3 = hi. Extending (f2"Tj) to the case 
cxp[A] cxp[_B] cxpL4] = cxp[Z] (see Section^}, it is found that K c g = 0, and thus 
it is retrieved that both bianisotropy and optical activity vanish in a medium 
with a center of symmetry [25[ . 

The present expansion in power series of frequency provides a new explana- 
tion for artificial magnetism and magnetoelectric coupling. Analytic expressions 
(|M|) of effective parameters can be used to analyze artificial properties. In parti- 
cular, we show from (1341) that : Artificial magnetism, previously proposed with 
high contrast [l|, [2j, llOj , can be obtained with arbitrarily low contrast ; and bia- 



nisotropy, previously achieved in J7-composites [27], can be present in simple 
one-dimensional multilayers. Note that, one can obtain more accurate asymp- 
totic expressions for the effective parameters with more terms in (|27|l and (IM)) , 
by taking higher order approximation in (1211) . 

Although this homogenized system has been studied by [28j, these authors 
assumed some magnetism and bianisotropy for the periodic multilayered stack, 
whereas in our case bianisotropy and magnetism come from a homogenization 
process (one might say ex nihilo). Moreover these authors assumed that the 
bianisotropy matrix was diagonal, which is not the case in the present paper. 
It is to the best of our knowledge the first time these constitutive relations 
are derived, and we emphasize that the mathematical theorem invoked in this 
section (BCH formula) can be used to generalize our result to two dimensional 
and three dimensional periodic structures [29j, such as woodpiles [llj. In the 
sequel, we shall also investigate numerically the stop band properties of such 
a periodic stack of dielectrics, and draw some illuminating parallels with the 



seminal paper by Pendry 12 1. 



4. Extension of BCH formula for m layers 

In Section [3l we have introduced our HOH algorithm for a periodic multi- 
layered stack consisting of an alternation of two dielectric layers, wherein the 
BCH formula is implemented. In this section, we would like to investigate the 
extension of HOH to a stack with m layers in a unit cell, correspondingly, a new 
form of BCH formula should be explored. We start with m = 3, which means 
a multilayered stack with an alternation of three dielectric layers is considered, 
the thickness of each layer in a unit cell is h m with m = 1, 2, 3, then the transfer 



10 



matrix of one unit cell will be 

T = T^Tj, = exp[iMifci] exp[iAf 2 /i 2 ] exp[zM 3 /i 3 ] (40) 

rewritten in a more general form, e.g. 

cxp[Z] = cxpL4i] exp[A 2 ] exp[A 3 ] (41) 

defines a product of three exponential functions. Obviously, it can be solved by 
an iteration of BCH formula. First, we suppose 

exp[A] = exp[Ai] expL4 2 ] (42) 

and A can be derived through equation (|2"Tj) 

A = A<® + A™ + A& + A^ + ■■■ (43) 

The term A^> represents the i th order approximation, and 

A(°) = A x + A 2 , 



A {2) = i^i, \A U A 2 \\ - -L[A 2 , IA U A 2 \ 



(44) 



then (|4Tj) turns to be 

exp[Z] = exp[A] exp[A 3 ] . (45) 

Using the BCH formula for Z : 



Z = log (cxpL4] exp[A 3 ]) = A + A 3 + -{A, As] 

+±[A, {A,A 3 jj ±{A 3 , [A,A 3 }] - 1[^3, {A, \A,A 3 \ 



(46) 



we suppose Z = Z'"' + Z^ + Z^> + ■ • • , where Z^ m ' is the m th order approxi- 
mation for Z. The zeroth order Z^ ' is simply the sum of A\, A 2 and A 3 , 

Z (o) = A (o) + A 3 = A! + A 2 + A 3 . (47) 

The first order including single commutator of these three matrices A\, A 2 and 
A 3 is 

Z {1) = A™ + \\AV>,M = \lAi,A 2 \ + l\Ax + A 2 , A 3 j . (48) 

and the second order including double commutators is 

Z {2) = AM + \\A^\A,\ + 1[A«, [AW.^l] - ±{A 3 , [A«, A 3 ]] 
= i[A ls [A 1( A 2 ]] - i[^ a , [Al A*}] + jiA 1; A 2 ], A 3 ] 



-L{Ai+A 2 , lA 1+ A 2 ,A 3 jj - l[As, [Ax + A 2) A 3 ]] . 



11 



(49) 



A similar algorithm holds for the third and higher orders, which will not be 
further explored here. 

Since the BCH formula for (|4"Tj) has been derived, one can easily realize the 
homogenization for a multilayered stack with an alternation of three dielectric 
layers. Here we assume that the third layer of the unit cell is identical to the 
first layer, i.e. £3 = e\ and /13 = hi, as well as M3 = Mi ; taking equation (|48[l 
with Ai = A3, we have Z^-> = 0. It should be noted that all the odd orders 
of approximation in the HOH asymptotics vanish, which can be attributed to 
the center symmetric property of the structure [30]. In contrast, even orders 
rule the approximation process in that case. Applying the formulae (|4"T1) - (|49I) 
to P0|) . one deduces the expressions for the effective parameters at 2nd order 
approximation 

£|l = 2ei/i + £2/2 — fJ>ofifa{ei ~ £2)(ei/i + £2/2) I 1 o ~ " 

3 \ uj z ^ £i£2 

uj 2 h 2 
e x = 2£i/i + £2/2 5—^0/1/2(21 - £2)(ei/i + £2/2) 

uj 2 h 2 
M|| = Mo + — — Mo/1/2 (ei - £2)(/i + fa) 

, ^h 2 2 / fc 2 £i +£ 2 
H± = fio + — 5— MoA/2(ei - £2)(/i + /2J 1 o 

S3- 1 = 2erVi + ejV2 + ^Mo/i/ 2 (£i - £2)(£rVi + ejVa) f 1 - ^J^l±^ 

3 V w Mo£i£2 

w 2 ft, 2 
Ms" 1 = Mo" 1 3— /i/2(£i - £2)(/i + fa) 

K±_ = K\\ = . 

(50) 
The effective bianisotropy K e s is equal to zero since Z^ 2p+1 ^ = 0, and only the 
artificial magnetism and high order corrections to the permittivity persist. Si- 
milar calculation can be applied to higher order approximation, e.g. the expres- 
sions of these parameters in 4th order approximation under a normal incidence 
is discussed in [3l|. 

So far, we have discussed the HOH asymptotic for a multilayered stack 
consisting of an alternation of two layers, as well as three layers ; and the BCH 
formula has been also amended correspondingly. If we extend this asymptotic 
procedure to a more general case, i.e. we consider a stack with an alternation 
of m(> 3) layers, then the transfer matrix of a unit cell becomes 

m 

T = exp[Z] = '[[exp[A i ]. (51) 

i=l 

Once again, tedious iteration of BCH in (fSTj) can produce all the formulae for 
different orders of approximation. Here, we just list the formulae from zeroth 
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(52) 



order to second order approximation : 

m 

z^ = E Ai , 

1 m i 

i=2 i=2 
^ m i— 2 m 

z(2) = iEiE^' A -ii<E^ 

1 m / i— 1 i — 1 i — 1 

These formulae can be checked by taking m — 3 and then compare with equa- 
tions (l4"Tl) - (|4l?)) . Apart from an iteration of BCH formula, another method to 
obtain the approximation for Z would be to expand each exponential function 
by Taylor series, and collect the terms with same order, which will not be further 
discussed in this paper. 

5. Corrector for HOH asymptotics 

In this section, we would like to introduce the corrector for the asymptotic 
error in the HOH algorithm, where a structure with thickness constant with 
respect to the frequency or the wavelength should be considered. We start with 
a structure consisting of two dielectric layers as shown in figure (2fa), where the 
parameters are e±, s 2l and thicknesses hi, h 2 satisfying h\ + h% ~ A. In order 
to obtain a homogeneous effective medium for such a structure, a geometric 
reconstruction is implemented here by dividing the structure into n times smaller 
unit cells, figure [2fb) shows the new construction when n = 2, the thicknesses 
of two layers being hi/2 and hi/2, respectively. It is noted that the thickness 
of each layer in a unit cell is decreasing in proportion to an increasing n, as 
shown in figure [He). When n tends to a large enough constant, the thickness 
of the unit cell hi/n will be much smaller than the wavelength of light, and a 
homogeneous effective medium can be achieved as shown in figure [5Jd) . 

The transfer matrices of the periodic medium and the effective homogenized 
medium should satisfy : 

{ exp[iM2h2/n] exp[iMihi/ri\\ c*. exp[iM e gh], for large n. (53) 

The BCH formula is still central to solve this problem, hence we recall its sta- 
tement : 

exp [ A x ] exp [ A 2 ] = exp [Ai + A 2 + -[Ai,A 2 } + —[^1, [^1, Mi 

-llA 2 ,lA 1 ,A 2 ]} + --]. (54) 



Let us take the first order estimate in (|54j) . This leads to the following error 
estimate : 
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Figure 2: Schematic diagram of the homogcnization algorithm for a structure with constant 
thickness : (a) A structure consisting of two homogeneous dielectric layers of permittivities 
ei, £2 and thicknesses h±, h,2 ; (b)-(c) A reconstruction of (a) by dividing the unit cell into 
n times smaller unit cells where n = 2 in (b), and a large n in (c) ; (d) An effective medium 
with thickness h. 



Proposition 5.1. If A\ and A 2 in (54\ ) are bounded by \\A\\ /2, then 



with 



||{exp [Ai/n] exp [A 2 /n]} n - exp [(A) + L A /n}\\ < — , 

n z 



(A)=A 1 +A 2 , L A = IA U Ail/2, 



a = ML exp [3 \\A\\] exp [||A|| 2 



Demonstration. Let S n and T n be denned by 

S„ = exp [(A) jn + L A /n 2 ] , 
T n = exp [Ai/n] exp [Ai/n] . 

One can write 

n 
nn rpn \ ^ qp— 1 ( Q 1^ \ r p n ~ P 



p=l 



so that 



i^-^ii<Eii 5 »ii p " 1 ii 5 »- t »iiii t » 



(55) 



(56) 



(57) 
(58) 
(59) 



0=1 
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Straightforward upper bounds for S n and T n are : 



||S„||<exp[||A||/n]exp[p|| 2 /n : 
||T„|| < exp[\\A\\ /n] . 
Then developing exponential functions in (|57p as series, 

S n = 1 + ((A) jn + L A /n 2 ) + ((A) /n + L A /n 2 ) 2 /2 + • • ■ 
T n = [ 1 + Ai/n + A 2 /(2n 2 ) + •••][ 1 + A 2 /n + A 2 2 /(2n 2 ) + 

one has 



■ ^^Lexp[2||A||]exp[||A|| 2 



3n 3 
Substituting §U§ and dHSJ) into ([SI]) provides the results flgJD - flgSp . 



(60) 



(61) 

(62) 
D 



In periodic classical homogenization [y, l32j | , only the leading order term is 
kept in the asymptotic procedure, hence the corrector is of order 1/n. Here, we 
find that : 

||{exp[Ai/n]exp[A 2 /n]} n -exp[<A)]|| <-, (63) 



with 



(A)=A 1 + A 2 , 6=||A|| 2 exp[2|| J 4|| 



(64) 



We therefore emphasize that our iterative procedure amounts to keeping more 
and more terms in the asymptotic expansion of classical homogenization and 
thus improves the order of the corrector of classical homogenization. Similar 
ideas have been implemented in the high-frequency homogenization recently 
developed by Craster et al. [9| , however with no correctors being derived therein. 
It would be also interesting to see how randomness would affect our correctors : 
at order zero, the corrector is known to vary between ■nT 1 ' 2 and n~ 



Furthermore, using the same proof process, we can also obtain the second order 
estimate of the limits (l54l . 



Proposition 5.2. If A\ and A 2 in |33) are bounded by \\A\\ /2, then 



|{exp [A x /n] exp [A 2 /n]} n - exp [(A) + L A /n + R A /n 2 ] || < — . 



with 



(A)=A 1 + A 2 , L A = \A 1 ,A 2 


1/2, 


A = ^[Ai, [AuAtH - 1^2,1^,^1 


\\A\\ 4 

■^-^exp[3 \\A\\] exp 


'\\M\ 


exp 


'\\Af_ 



(65) 



(66) 
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Demonstration. Let S n and T n be denned by 

S n = exp [(A) jn + L A /n 2 + R A /n 3 ] 
T n = exp [Ai/n] exp [A2/71] . 

Upper bounds for S n and T n are straightforward : 

||5 n || < exp [||A|| /n] exp ||A|j 2 /n 2 exp 
||T„||<exp[||A||/n] . 
And developing the exponential function in (|67p as a series , 

S n = l + ((A) jn + L A /n 2 + R A /n 3 ) + ((A) /n + L A /n 2 + R A /n 3 ) 2 /2 + ■ 
T n = [ 1 + A 1 /n + A 2 /(2n 2 ) + •••][ 1 + A 2 /n + A 2 /(2n 2 ) + • • • ] 

one has 

Mil 4 



(67) 



(68) 



(69) 



O _ T II < 



4n 4 



exp[2||A||]exp \\A\ 



exp 



\A\\ 



Substituting (ggj, dTOJ) into ([59]) leads to (|65|) -(|66 



(70) 



□ 



It is noted that as the approximation order increases, the speed of the conver- 
gence defined by the difference between transfer matrices of multilayers and ef- 
fective medium increases by a factor 1/n, hence it seems natural to conjecture 
that for the higher order approximation, the estimate between the transfer ma- 
trices of multilayers and the effective medium will be much more accurate with 
an error of l/n p , with p the order taken in HOH approximation process. 

Similarly, the asymptotic corrector can be applied to the stack with m layers. 
Here, we explore the corrector for HOH approximation in a multilayered stack 
with three layers, where the permittivities are E\, e 2l £3 and thicknesses are 
hi/n, hz/n, h^/n. Applying the same geometric reconstruction as shown in 
figure [21 the equivalent relation between the transfer matrices of the stack and 
effective medium is 



{ exp[iMihi/n] exp[iAf2^-2/^] exp[iM3/i3/n]} ~ exp[iM e gh] . 

According to equations (gDJ) and ([37 ]) - P? ]) . the BCH leads to 

exp [A x ] exp[A 2 ] exp[A 3 ] = exp [A x + A 2 + A 3 + {A X A 2 - A 2 A 1 )/2 

+ {A X A Z - A 3 A x )/2 + {A 2 A 3 - A 3 A 2 )/2 + 

Taking the zeroth order approximation as an example, we can state 



(71) 



(72) 
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Proposition 5.3. If A\, A 2 and A 3 in J7g[ ) are bounded by \\A\\ /3, th 



en 



a 



(76) 



||{exp[A 1 /n]exp[A 2 /n]e X p[A 3 /n]} n -exp[<A)]|| < -, (73) 

n 

with 

(A) =A 1 +A 2 + A 3 , 

(74) 
a= ||A|| 2 exp[2 \\A\\] . 

Demonstration. Let S n and T n be denned by 

S n = cxp [{A) /n] 

(75) 
T„ = cxp [Ai/n] exp [A 2 /n] exp [A 3 /n] . 

Equations (|58p and (|59[) remain valid and upper bounds for S n and T n are 

||5 n ||<exp[||^||/n], 
||T„|| < exp [|| A|| /n] . 
Then developing the exponential function as a series, 

S n = 1 + ((A) jn) + ((A) /nf /2 + • • • 
T n = [ 1 + A x /n + A\/(2n 2 ) + •••][ 1 + A 2 /n + A 2 2 /(2n 2 ) + • • • ] (77) 

[l + A 3 / n + Al/(2n 2 ) + ---] , 

one has 

||S„-T n ||<MLexp[2||A||]. (78) 

Substituting (J7SJ) and (75) into (59) leads to (73) and (75). D 

The proof indicates that the corrector is in order of n^ 1 when taking the 
classical homogenization (zeroth order approximation) for a multilayered stack 
with an alternation of three layers, this can be adopted for the m layers case. 
Correctors for higher order approximation, as well as for a multilayered stack 
consisting of an alternation of m layers can be obtained by the same algorithm. 

6. Numerical calculations : Dispersion law and transmission curves 

In this section, we would like to numerically investigate the asymptotic de- 
gree between the multilayered stack and its effective medium obtained from 
HOH algorithm, where the dispersion law and transmission curves are explored. 
According to the previous analysis, the transfer matrix defined by the exponen- 
tial function of matrix M is analytic, and it can be expanded as a Taylor series, 
taking the effective transfer matrix as an example, 






T cff = expHMeff/o = ^(-i)*-^- + 2,Hr +1 ;° 1 -, • (79) 



n 



Considering a s-polarized incident wave, the column vector in ([7]) is denned by 



F=[E ± ,H|,] T 
The matrix M c g in (|15|) is a 2 by 2 matrix, and 



Ks = 



iuK 



-we j 



k z 



— W/i|| 

-icoK 



2 












1 


n 




— h 2 





i 








- 1 



where 



M3 



(80) 



(81) 



(82) 



Plugging (1811) into (J79I) and considering Taylor series of the sin — cos functions 



MA) = ^ J J „. A 2 " +] 



n=0 



(2n + l)!' 



cos 



M = £ 



n=0 



(~l) r 

(2n)l 



T 2 



(83) 



we obtain 



.M, 



off 



T eff = cos(k cS h) - i- sin(fc e ff/i) 

«eff 



cos(fc c ff /i) + wi^i 



sin(fc off /i) 



sin(fc eff /i) 



fceff 

A; 2 1 sin(fc e ff/i) 
j( W£j _ ) — 

W /i3 Keff 



IWjUll - 



COs(fc c ff/l) — w-RT. 



s'm(k cS h) 



ivir 



Similarly, the transfer matrix of the dielectric layer C m is 

sin(0 m h m ) 



(84) 



T„ 



cos (f3 m h m ) 
k 2 1 sin(/3 m ft m ) 



iW/iQ- 



A, 



i(we m ) 

w Mo 



A, 



cos ((3 m h m ) 



« 2 _ 2 i 2 



(85) 

The transfer matrix T of the unit cell consisting of two dielectric layers is derived 
from the above expression as 

T = T 2 Ti . (86) 

6.1. Dispersion law 

A general expression for the dispersion law in a periodic structure is defined 
by the trace of transfer matrix T of a single period [34j, |3a, [36j. Since the 
eigenvalues and eigenvectors of T (and thus any power of T) are the Bloch wave 
vectors and Bloch states of the periodic structure (in the limit of infinite n), 
further physical insight can be achieved in the single period matrix T . Hence, 
for a multilayered stack with two layers, we have 

tr(T)/2 = cos (fi!hi) cos (f3 2 h 2 ) - h£l + E±) sin {p x h x ) sin (fch 2 ) (87) 

£ Pi Pi 
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while for the effective medium, 



tr(T eff )/2 = cos(k eS h) . 



(88) 



Here, k e g defined by l|82p can be obtained by substituting the expressions of the 
effective permittivity, permeability and bianisotropy, which are derived from the 
HOH algorithm. 

Considering a normal incident plane wave (fc = 0), the s-polarization and 
p-polarization coincide, since e» = £j_, /in = /Uj_, Ku = K±. Hence we take a 
s-polarized incident wave as an example, and assume the two dielectric layers 
of the stack are Glass and Silicon, respectively ; the relative permittivities are 
£1 = 2, e 2 = 12 and the filling fraction are /i = 0.8, fa = 0.2. For the sake 
of illustration, in the HOH algorithm, we take the 3rd, 7th and 19th order 
approximation for the effective medium. The curves of the effective permitti- 
vity, permeability and bianisotropy in 19th order approximation versus norma- 
lized frequency are depicted in figure Eta), where u) — ujh/c, the expressions of 
these effective parameters are omitted here to save space. It is observed that all 
three curves are increasing along with the frequency, wherein the effective per- 
meability (dash red line) has values greater than 1, and effective bianisotropy 
(dotted-dash blue line) is non- vanishing. In other words, artificial magnetism 
and bianisotropy can be achieved from dielectrics through HOH asymptotics, 
as it has been predicted in the theoretical analysis of Section [31 
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Figure 3: (a) The curves of effective permittivity (solid black line), permeability (dash red 
line) and bianisotropy (dotted-dash blue line) in 19th order approximation ; (b) Dispersion 
laws of the multilayered stack (solid black line) and its effective medium in 3rd order (dotted 
red line), 7th order (dash blue line) and 19th order (dotted-dash black line). HOH breaks 



down at i 



: 0.202, lower edge of the first stop band. 



Figure GUb) shows the dispersion law of the stack, as well as that of the 
effective medium in different order approximations, which is obtained by sub- 
stituting the expressions of the effective parameters into k e s of (f8"2"]> and then 
((55)) . The lower edge of the first stop band is denoted by & g = 0.202, where 
tr(T)/2 = — 1 [37]. It should be noted that a good agreement between the dis- 
persion laws of the stack and effective medium can be observed at the lower 
frequency band, and the asymptotics of these two curves can be improved by 
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taking higher orders of approximation, e.g. 19th order approximation (dotted- 
dash line) shows a better asymptotic with the multilayers (solid line) from zero 
frequency (quasi-static limit) up to a normalized frequency around 0.18. 

Moreover, an oblique incident plane wave in s-polarization as well as in p- 
polarization are also analyzed, wherein the same parameters of the stack are 
taken and the incident angle is 8i = 30°, the dispersion laws of the multilayered 
stack and the effective medium are shown in figure 2J The 3rd and 7th order 
approximations are considered in the HOH algorithm, similarly, asymptotics 
improve in conjunction with higher orders of approximation for both s- and 
p-polarized waves. 





(a) s-polarization 



(b) p-polarization 



Figure 4: Dispersion laws of the stack (solid black line) and its effective medium in 3rd order 
(dotted red line), and 7th order (dash blue line) approximation, under an oblique incident 
wave with 0i = 30°, (a) s-polarization with Ci g = 0.206, (b) p-polarization with u) g = 0.215. 



6.2. Transmission curve 

Apart from the dispersion law, asymptotics between the transmission curves 
of the multilayered stack and the effective medium is another important fea- 
ture to be checked. Figure [5] shows a schematic diagram of the reflection and 
transmission for an incident wave U z on the effective medium, the thickness of 
which is denoted by nh, the upper and lower spaces of the effective medium are 
supposed to be vacuum with permittivity £n and permeability [Aq. 

We assume the incident wave in form of Fourier decomposition in ([3j) is 



U =U exp[-ife 3 a;3], x 3 > 0, 



(89) 



where fc| = w 2 £oA t o ~ kf — k%, and Uo the amplitude of the incident electroma- 
gnetic field. The reflection and transmission waves are 



U = U rexp[jfc 3 a;3], 



x 3 >0 



U = Uo £exp[— ika(x3 +nh)], x 3 < —nh 



(90) 



20 



L" 

\ 



nh 



/ 



\. 



u' 



Figure 5: Schematic diagram of the reflection and transmission for an incident wave on a 
periodic media with thickness nh, the upper and lower spaces are vacuum. 



For a polarizable incident wave, the column vectors F defined in ([7]) are 

s — polarization : F = [Ej_,H|| ] T 
p — polarization : F — [ En , H^ ] T 



F = 



(91) 



while for the upper and lower vacuum of the effective medium, (|9ip will be 
simplified as 

U 



with 



(92) 
(93) 



s — polarization : U = Ej_ , vq = uq , 
p — polarization : U = Hj^ , vq = Sq . 

Assume T n h is the transfer matrix of the effective medium with thickness nh, 
we have 

F(0) = T nh F(-nh) (94) 

with 



Tn 



rpTi 



h\C n -i{a) - C n - 2 {a) t 12 C n -i{a) 
t2iC n -i(a) tmCn-iia) - (7 n _ 2 (a) 



(95) 



where a = (in +t22)/2 = tr(T)/2, T the transfer matrix of unit cell. And C n is 
the Chebyshev polynominals of the second kind 38] 



C n {a) 



sin[(n + l)cos 1 a] 



(96) 



Again, if we consider a s-polarized normal incidence, then U = E^ and H|| 
i(ujuo)~ 1 d X:i E±, (|94|) turns to be 



1 + r 

(1-r) 

W/XQ 



= T. 



(s) 



t 

UJUo 



(97) 
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and the transmission coefficient is 



(r« 



■T& ) )/2+[T® + 



(s) 



(fc3/( W/ io)) 2 T^]/2 



(98) 



According to the duality between s- and p- 



-.(p) 



replace ^o by eo, as well as T^ by T^' in 
coefficient for p-polarization, e.g. 



polarization, one only needs to 
(l98l) to obtain the transmission 



h = 



(T&> + Tfr>)/2 + [T. 



(p) 
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(VM) 2 ^]^ 



(99) 



Let us consider again a multilayered stack with e\ — 2, e 2 = 12, /i = 0.8 and 
/a = 0.2 as an example, the thickness of the structure is supposed to be nh 
with n a constant, i.e. n = 20, to allow a numerical calculation in Matlab. Note 
that the s- and p-polarized incident waves coincide under a normal incidence, 
i.e. t s = t p = t. Applying T = T c g as shown in (|84[) . and T = T-£T\ with T m 
in (|85p to (j9"5)) and (|98p , the transmission curves of the effective medium in 7th 
(dotted-dash red line) and 19th (dash blue line) order approximation are shown 
in figure El as well as the transmission curve of the stack (solid line). 

It can be observed that the lower order approximation (dotted-dash red 
line) only fits well with the curve of the multilayer (solid black line) in the 
low frequency band ; while an improved asymptotic (dash blue line) can be 
achieved by higher order approximation (e.g. 19th order). In agreement with 
the dispersion law, the asymptotics between the two transmission curves of 
multilayer (solid line) and effective medium in 19th order approximation (dash 
line) become invalid near the lower edge of the first stop band. 




0.1 0.15 

a/(2ir)(«M 



Figure 6: Transmission curves of the multilayered stack (solid line) and effective medium in 
7th order (dotted-dash red line) and 19th order (dash blue line) of approximation. 

Similarly, the transmission curves of the stack and effective medium under 
an oblique incidence with 9i = 30° in s-polarization, as well as p-polarization are 
shown in figure [7] Once more, the asymptotic approximation between the two 
transmission curves of the stack and the effective medium is quite good in the 
low frequency band, and improves with higher order approximation in effective 
medium as shown in the dispersion laws of figure |4] 
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(b) p-polarization 



Figure 7: Transmission curves of the multilayered stack (solid line) and effective medium 
in 3rd order (dotted-dash red line) and 7th order (dash blue line) of approximation : (a) 
s-polarization, (b) p-polarization 



6. 3. On logarithm of transfer matrix and analyticity 

From figures [4][7l it should be noted that asymptotics break down around 
the lower edge of the first stop band dig between the dispersion laws and the 
transmission curves of the multilayer and its effective medium, no matter how 
high the order of approximation is. This invalidity is contributed to the po- 
wer series expansion of X = iM e f[h in (|2ip which diverges at (2> g . Indeed, we 
choose BCH formula to obtain the approximation for matrix M c g furthermore 
for effective permittivity, permeability and bianisotropy, where we have taken 
X = log{exp[A] exp[i?]} in ([2Tjl . In complex analysis, a branch of log(z) is a 
continuous function L{z) defined on a connected open subset G of the complex 
plane, such that L(z) is a logarithm of z for each z in G [39|. An open subset G 
is chosen as the set C — M<o obtained by removing the branch cut (thick solid 
line) along the negative real axis and the branch point (empty point) z = from 
the complex plane, as shown in figure E^b). 

On the other hand, the transfer matrix T — exp[A] exp[B] can be factorized 
by eigendecomposition as 

T = QAQ- 1 (100) 

where Q is a square matrix consisting of the eigenvectors of T, and A is a 
diagonal matrix whose components are the eigenvalues (denoted by A±), and 



X, 



a ± i\f 1 — a 2 



(101) 



with a = tr(T)/2 the half trace of transfer matrix T . Furthermore, in order 
to derive the effective matrix M e g for those effective parameters, we take the 
logarithm of T 



AL 



log{T} = Qdiag[log(A + ),log(A_)]Q- 1 . 



(102) 



The curve of tr(T)/2 versus frequency is depicted in figure (5Ja), and one has 
A + = A_ = a = +1 at the zero frequency (denoted by cross sign), while 
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Re(z) 



Figure 8: (a) Curve of tr(T)/2 versus frequency for the mutlilayered stack consisting of an 
alternation of two dielectric layers ; the eigenvalues of the transfer matrix are denoted by \± 
and vary from +1 (cross sign) to —1 (filled point) in the lower pass band ; while tr(T)/2 = —1 
defines the edge of the first stop band ; (b) Schematic diagram of a branch of log(.z), a region 
G (connected open subset of the complex plane) can be typically obtained by removing from 
C the interval (— oo,0]. For the logarithm of X±, we choose the upper and lower half circle 
paths for \± to approach —1 at first edge of stop band (filled point) lying on the negative 
real axis. 



A + = A_ = a = — 1 at the edge of the first stop band Cj g (denoted by filled 
point). In the lower frequency band, we choose A + and A_ (starting from +1) 
which approach — 1 by the upper half-circle path and the lower half-circle path, 
respectively, where the paths lie in the open subset G : log(A±) is always analy- 
tic and unique. However, at the edge of the stop band where A + = A_=a = — 1 
on the negative real axis, the logarithm is no longer analytic when its arguments 
meet at the branch cut of the logarithm : This implies that an expression of M e g 
as a power series of the frequency u has its radius of convergence bounded by u> g . 
In other words, the effective parameters lose their efficiency for the asymptotic 
approximation at frequencies higher than the first stop band, but they work just 
fine in the lower pass band. In order to achieve all frequency homogenization 
for a periodic structure, a new set of effective parameters (e.g. effective refrac- 
tive index and surface impedance) should be introduced [3l| , where the analytic 
property of transfer matrix in the complex plane is ensured. 



7. Frequency power expansion of the transfer matrix 

Although the function log{exp[A] exp[B]} is no longer analytic at the lower 
edge of the stop band, the transfer matrix T = exp[A] exp[_B] is analytic in the 
whole complex plane and can be approached by a power series at any frequency, 
a fact which will be numerically checked in this section. In mathematics, an 
exponential function can be approximated by a Taylor series as 



A n 



ex P^] = ^17 = ! 



A 



n=0 



A*_ 

2! 



3! 



(103) 



24 



Hence, the transfer matrix T takes the form 
T = T X T 2 = exp[iMihi] exp[iM 2 h 2 ] = exp[wiVi] exp[u}N 2 ] 

UmV , (^^i) 2 , W) 3 , \[, , AT , Mfe) 2 , ( W iV 2 ) 3 

1 + uJVi H — 1 — 1 1 + ujN 2 + 



2! 3! / V 2! 3! 

(104) 
with notation wiV, = iMihi. Let us expand and organize (|104|) in powers of w, 

/V 2 /V 2 
T c ff = TiT 2 = 1 + w(JVi + N 2 ) + w 2 (^± + -J. + NlN2 ) + ... (105) 

We collect the terms from ui° to uj p in (|105|) as an ansatz for T e g 

T cff ~ T (0) + wT (1) + w 2 T (2) + w 3 T (3) + ■ ■ • + w v T {p) . (106) 

Obviously, with increasing p, the approximation in (|106[) becomes more accurate. 

Considering a normal incident wave in s-polarization, the matrices N m read 

as 

/i 
-£m 



Nm. = ih r 



1, 2 (107) 



so that substituting them into (|106|l the half trace of T e ff can be expressed as 

~ 2 ~ 4 "4 

tr(T cff )/2 = 1 - y (ei/i + £2/2) + ^(£1/1 + e 2 / 2 ) 2 - ^(ei - E a ) a /?/l 

-^(£1/1 + £ 2 / 2 ) 3 + ^/i 2 /!(£i - £ 2 ) 2 (/i - / 2 )(ei/i - £2/2) 

+ ^/i 2 / 2 2 (£i ~ e 2 ) 2 (ei/i + £2/2) + ^A 3 /|(ei - e 2 ) 2 (ei + £2) 

+^/i /l(ei " £2) [(el/a - £ 2 /i) + (£2/1 - e 1 f 2 )(e 1 f 1 + e 2 f 2 )] + • ■ 

(108) 
where the normalized frequency u = u)h/c. It is noted that there are only terms 
containing even order of ui, which is due to the fact that the T( 2 p +1 ) in (|106| l 
are all off-diagonal matrices with zero diagonal components. 

The curves of tr(T)/2 versus frequency are shown in figure [Sf a) : The thick 
solid line represents the half trace of transfer matrix for a multilayer consisting 
of an alternation of two dielectric layers with same parameters as assumed in 
Section [51 the dash green line, dotted red line, dotted-dash blue line and the 
thin solid line are tr(T e g)/2 with p = 2, 4, 6, 8 in (|106l) for T c s , respectively. 

By comparison with the curve of the multilayered stack (the thick solid line), 
it can be observed that the approximation for T e g in (|106p with p — 2 (dash 
green line) is just efficient in a range of low frequency, while p — 4 (dotted line) 
gives a sharper estimate for the multilayered stack, but it totally misses the stop 
band. Moreover, if we push the approximation to p = 6, the dispersion curves 
are nearly superimposed up to the edge of the first stop band, so well beyond 
the range of validity of classical homogenization [6] . However, the approximation 
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0.2 0.3 

ffl/(2jr)(«,=0) 



(a) 



(b) 



Figure 9: (a) Curves of tr(T)/2 versus frequency for the multilayered stack (thick solid line) 
and the effective medium under a normal s-polarization incidence, the dotted, dash, dotted- 
dash and thin solid lines correspond to the approximation for T c g in I I1U6H by taking p = 
2,4,6,8. (b) Transmission curves of the multilayered stack (solid line) and the effective medium 
(dash line) with p = 8 in ||106|I . Parameters of the dielectric layers of the stack are assumed 
to be the same as in Section [6] 



with p — 6 (dotted-dash curve, which is always decreasing) breaks down at the 
lower edge of the stop band. In order to better approximate T e g- , one needs to 
push the approximation to the next even number of p (i.e. p = 8, the thin solid 
line), which changes the curvature and gives a sharper estimate in the stop band 
region, although its intersection with the horizontal axis defines an approximate 
position for the upper edge of the stop band. This can be improved by taking 
larger p in (|106|) . Altogether, the larger p taken in (|106[) . the more accurate the 
approximation between the dispersion curves of the effective medium and that 
of the multilayered stack. 

Moreover, according to the expression of the transmission coefficient in (|98|) . 
we take p — 8 in (|106[) for T e g, the transmission curves for both the multilayer 
(solid line) and the effective medium (dash line) are depicted in figure Mjo). 
A good agreement between these two curves can be observed up to the first 
stop band, as predicted in figure [SJa). Similar calculation can be applied to an 
oblique incidence. This demonstrates that the transfer matrix of the effective 
medium can be approached as a frequency power series at any frequency. 



8. Concluding remarks 

We provide a rigorous high-order homogenization (HOH) algorithm for one- 
dimensional moderate contrast photonic crystals, where the period of the struc- 
ture approaches the wavelength of optical waves. From an expression of transfer 
matrices in terms of exponential functions, S. Lie and BCH formulae are applied 
in the HOH asymptotic. Analytic expressions of the effective parameters are de- 
rived for a stack with two layers in Section [3J where the artificial magnetism 
and magnetoelectric coupling effect are achieved in such a moderate contrast 
periodic structure. In Section^ we explore the extension of HOH algorithm to a 
stack with an alternation of m dielectric layers, and derive the expressions of the 
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effective parameters for a center symmetric stack : The magnetoelectric coupling 
vanishes while the artificial magnetism can be achieved with non-resonant per- 
iodic structures. Furthermore, the corrector for the asymptotic approximation 
of a finite stack by its effective medium has been discussed in Section [5] : The 
asymptotic error is of order l/n p with p the order of the approximation. Finally, 
based on the expressions of the effective parameters, we numerically validate our 
approximation method by comparing both the dispersion law and transmission 
property of the stack and its effective medium in Section |H1 The good agree- 
ment between these curves demonstrates that the asymptotic approximation is 
efficient throughout the lower pass band, while at the edge of first stop band 
the logarithm function is no longer analytic. Finally, we investigate the approxi- 
mation for the transfer matrix instead of matrix M c g of the effective medium 
by a frequency power expansions, the dispersion law as well as the transmission 
curves of the transfer matrices for both the multialyer and the effective medium 
are explored in Section[71 and the excellent agreement confirms that the effective 
transfer matrix can be approached by a power series at any frequency. 
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